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Abstract. Massive disk fragmentation has been suggested to be one of the mechanisms leading to the formation of giant 
planets. While it has been heavily studied in quiescent hydrodynamic disks, the effect of MHD turbulence arising from the 
magnetorotational instability (MRI) has never been investigated. This paper fills this gap and presents 3D numerical simulations 
of the evolution of locally isothermal, massive and magnetized disks. In the absence of magnetic fields, a laminar disk fragments 
and clumps are formed due to the effect of self-gravity. Although they disapear in less than a dynamical timescale in the 
simulations because of the limited numerical resolution, various diagnostics suggest that they should survive and form giant 
planets in real disks. When the disk is magnetized, it becomes turbulent at the same time as gravitational instabilities develop. 
At intermediate resolution, no fragmentation is observed in these turbulent models, while a large number of fragments appear 
in the equivalent hydrodynamical runs. This is because MHD turbulence reduces the strength of the gravitational instability. As 
the resolution is increased, the most unstable wavelengths of the MRI are better resolved and small scale angular momentum 
transport starts to play a role: fragments are found to form in massive and turbulent disks in that case. All of these results 
indicate that there is a complicated interaction between gravitational instabilities and MHD turbulence that influences disk 
fragmentation processes. 
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1. Introduction 

The recent discovery of more than 100 extrasolar planets has 
renewed interest in planet formation theory. Because of se- 
lection effects, all of the planet discovered so far are planets 
with masses larger than the mass of Neptune. However, the 
problem of their formation remains unsolved, despite the at- 
tention it continues to enjoy in the r ecent literature jBossll998l 
iMaver et alJ2002tlRice et alJ2003l) . Two very different models 
have been proposed. 

In the so-called "core accretion model," giant planets first 
form their cores by dust accumulation. When the core reaches 
a mass of the order of 10 earth masses, a phase of rapid 
gas accretion ensue s and builds up the envelope of the planet 
(Poll ack etaill99d) . The typical timescale to form a planet of 
1 Jupiter mass in a realistic disk is of the order of 10 7 years. 
This rather long timescale has raised questions since it is simi- 
lar to the expected lifetime of the disk. However, recent studies 
suggest that the migrat ion of the planet in the disk may over- 
come this difficulty jHourigan & Wardll984HRice & Arm itage 
l2003UAlibertetalJ2004 . 



The alternative scenario is that giant planets form directly 
in massive protopl anetary disks by gravitational fragmentation 
llBossll 19981 EoO O*). Massive disks are linearly unstable to ax- 
isymmetr ic gravitation al instabilities when the Toomre Q pa- 
rameter ( Toomr eTl964l) . defined by 



Q 



ttGZ 



(1) 
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becomes smaller than 1. In the definition above, c s is 
the sound speed, k is t he epicyclic frequency (see, e.g., 
iBinnev & Trernainelll987l) . 2 is the disk surface density and 
G is the gravitational constant. Gaseous disks are also unsta- 
ble to non-axisymmetric perturbations when Q > 1 . In the later 
case, spiral arms develop and transport angular momentum out- 
ward in the disk. If the nonlinear evolution of the instability 
is very violent, these spiral arms may fragment to form dense 
clumps which have been in terpreted as proto-giant planets. For 
example, ICiammiel (hoOl) found that a disk fragments when 
its cooling ti mescale is shorter than t coo \ ~ 3Q~'. This was 
confirmed by Rice et al. (2003) using SPH numerical simula- 
tions. lMaver et alJ j2002l) : lMaver et all J2004 also found frag- 
mentation fo£_avarietv of dis k models and equations of state. 
However, Pic kett et alJ (2000) failed to obtain fragmentation 
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in similar thermodynamic conditions. Whether or not massive 
disks fragment remains controversial. 

As described above, the issue of disk fragmentation has 
been largely addressed in the past few years. These studies omit 
the effect of magnetic forces, despite the known effect of even 
a weak magnetic field on the dynamical evolution of a disk. 
Indeed, rotating disks in which the angula r velocity decreases 
outward are unstable to the MRI jBalbus & Hawlevlll99lh . In 
the last decade, numerical simulations, both local and global, 
have shown that its nonlinear outcome is MHD turbulence that 
transp ort s angular mom entum outward (see Balbus & Hqwlg^ 
1998, or Balbus 2003, for a review). Recently, Fromang et alJ 
(2004a, 2004b) investigated the evolution of massive, magne- 
tized disks, in which the MRI and gravitational instabilities are 
expected to develop simultaneously. The main result of this 
study is that both instabilities are tighly coupled and strongly 
interact with one another, leading to a quasi-periodic variation 
of the mass accretion rate onto the central star. However, these 
results were obtained using an adiabatic equation of state, a 
regime in which no gravitational collapse was found in the disk. 
The latter is more likely to occur in isothermal disks, since the 
pressure support is less important in that case. The main goal 
of this paper is to investigate the influence of MHD turbulence 
on the evolution of isothermal disks, with an emphasis on the 
question of their fragmentation. 

The plan of the paper is as follows. Section 2 presents 
the numerical methods used and the initial conditions of the 
simulations. Section 3 describes their results, focusing first on 
purely hydrodynamic cases and then on full MHD simulations. 
These results are finally discussed in section 4. 

2. Method 

The relevant equations are those of ideal MHD for a self- 
gravitating fluid, in which p is the mass density, v is the ve- 
locity, B is the magnetic field, P is the gas pressure and O = 
<$>s + $>c is the total gravitational potential, which has contri- 
butions <t>s from the disk self-gravity and <t> c from a central 
mass: 



is adiabatic, the pressure is derived from the internal energy e 
using 



% + V-(pv) = 0, 

of 



(2) 



tdy \ 1 

p — + (T-V)v = -VP - pVO + — (VxB)xB, (3) 

\ at I 4k 



^ +V - V )S= VX(VXB) ' 

V 2 «D S = 4nGp, 



(4) 
(5) 



The equation of state is either adiabatic or locally isothermal. 
In the latter case, the temperature in the disk is a function of 
position, but not time: 



- =f(r,(f>,z). 
P 



(6) 



P = (y-l)e 



(7) 



where y — 5/3 is the poly tropic index, e is calculated using the 
following energy equation: 



(8) 



P \8t + ^)\-p) = - P ^ 



Here (r, <p, z) denotes the standard cylindrical coordinates that 
will be used throughout the paper. When the equation of state 



2.1. Numerical procedure 

To solve t he equations described above, I used the code 
GLOBAL jHawlev & Stond 1 1995b . It uses time-explicit 
Eulerian finite differences and is designed to handle stan- 
dard cylindrical coordinates. The magnetic field is evolved us- 
ing the combined Method of Characteristics and Constrained 
Transport algorithm (MOC-CT), which preserves the diver- 
gence of the magnetic field to machine accuracy at all time. 
In order to reduce the time of the simulations, the computa- 
tional domainJ^_jeduced_to [0, 7r] in the azimuthal direction, 
as in iFromang et alJ d2004bl) . Finally, outflow boundary con- 
ditions are applied in the radial and vertical directions, while 
periodic boundary conditions are applied in (f>. 

The self-gravitating potential O s is calculated with 
a Poisson solver we ll adapted to cylindrical geometry 
(Froman g et alJl2004b) . The density distribution of the disk is 
first Fourier transformed in the azimuthal direction. For com- 
putational reasons, only those Fourier components p m whose 
index m ranges between and m max are retained at this stage. 
For each of them, the procedure is the following: (1) From 
p m , the gravitational potential is obtained at the boundary 
of the computatio nal domain usi ng an expansion in Legendre 
functions JCohl & Tohlrnel fr999). (2) is calculated every- 
where on the comp utational doma in using the Successive Over 
Relaxation method (Hirsch 1988). Finally, O s is reconstructed 
using the values of O™, with m e [0, m max ]. 

One has to understand that the value of m max described 
above is crucial for the problem of disk fragmentation ad- 
dressed in this paper. Indeed, the larger the number of coef- 
ficient included in the calculation of <P$, the smaller the struc- 
tures that can be described by the simulation. Of course, this 
comes also with a larger computational time. Nevertheless, in 
order to accurately describe the tiny structures associated with 
the possible fragmentation of the disks, I used m max = 32, 64 
or 128 in this wor k, in contrast with the lower value m max = 8 
typically used by Froman g et alJ d2004bl b) in adiabatic simu- 
lations where only large structures appeared. 

2.2. Initial conditions 

The properties of the various runs performed are summarized 
in table [2 The first column gives the label of the models. HD 
refers to hydrodynamic models and T refers to runs starting 
with an initial toroidal field. Column 2 gives the number m max 
of Fourier coefficient included in the calculation of the grav- 
itational potential. The third column gives the initial ratio (f3) 
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Model 




<fi> 


Resolution 


EQS 


Fragments? 


HDIS02 


32 


+ oo 


(128,64, 128) 


iso 


No 


HDIS03 


64 


+ 0O 


(128, 128, 128) 


iso 


Yes 


HDIS04 


128 


+ 0O 


(128,256, 64) 


iso 


Yes 


T2* 





8 


(128, 64, 128) 


adia 


No 


T3* 





8 


(128, 128, 128) 


adia 


No 


T4* 





8 


(128, 256, 64) 


adia 


No 


TIS02 


32 


8 


(128,64, 128) 


iso 


No 


TIS03 


64 


8 


(128, 128, 128) 


iso 


No 


TIS04 


128 


8 


(128,256, 64) 


iso 


Yes 



Table 1. The first column gives the model label. The second 
and third columns respectively give the number m ma x of Fourier 
coefficient included in the calculation of the gravitational po- 
tential and (J3), the initial ratio of thermal to magnetic pressure. 
The resolution (N r> Nf,N z ) of the model is shown in column 4, 
while column 5 gives the equation of state used ("iso" stands 
for locally isothermal and "adia" for an adiabatic equation of 
state). Finally, the last column indicates whether or not frag- 
ments were found in the isothermal runs. 



of the volume averaged thermal pressure to the volume aver- 
aged magnetic pressure. Column 4 gives the number of grid 
points (N r , N,/,, N z ) of the run and column 5 describes the equa- 
tion of state ("iso" refers to a locally isothermal equation of 
state). Whether or not fragments form in the disk is described 
in the last column. 

Some care has to be taken regarding the initial conditions 
of each of these models . Model T2* was already presented 
in Froman g et alJ rt2004bh : assuming an adiabatic equation of 
state, it calculates the evolution of a disk initially in equilib- 
rium, with a mass half that of the central mass, threaded by 
a toroidal weak magnetic field for which the ratio of the vol- 
ume averaged thermal pressure to the volume averaged mag- 
netic pressure (J3) is initially 8. The MRJ grows because of the 
magnetic field, and the disk eventually becomes fully turbu- 
lent. No gravitational instability develops since m max — in 
this model, which is evolved for about 6 orbits at the loca- 
tion of the initial outer edge of the disk. At this stage, it has 
reached a quasi-stationary turbulent state typical of the out- 
come of global simulations of zero mass disks: mass slowly 
falls onto the central object under the action of the Maxwell 
and the Reynolds stresses, while angular momentum is trans- 
ported in the outer parts of the disks, which spreads as a conse- 
quence. At the same time, a tenuous, magnetized corona builds 
up above and below the disk main body. 

This quasi-stationary state emerging from model T2* after 
6 orbits was used as the initial condition for model TIS02. To 
do so, the temperature was calculated at the end of model T2* 
everywhere on the grid and then kept constant for the further 
evolution of model TIS02, whose equation of state is locally 
isothermal. 

In order to compare model TIS02 with an hydrodynamic 
model whose properties are as closely matched as possible, the 
disk structure resulting from model T2* was also used to com- 
pute the initial conditions for model HDIS02. The procedure 
used to do so is as follows: the result of model T2* was first 



Fig. 1. Snapshots of the logarithm of the equatorial den- 
sity of the disk for model HD2, respectively at times t = 
9.0, 9.2, 9.3, 9.5, 9.9 and 10.3 from top left to bottom right. The 
structure of the disk varies very rapidely. Short lived high den- 
sity clumps are visible in each snapshots. 



azimuthally averaged and the magnetic fields were removed. A 
small random perturbation was then applied to the disk before 
restarting the simulation with a locally isothermal equation of 
state. Note that this sudden removal of the relatively weak mag- 
netic field does not strongly affect the disk global structure. 

The other simulations were initialized in a similar way. 
Model T3\ TIS03, HDIS03 and T4*, TIS04, HDIS04 are 
the high resolution equivalent to model T2*, TIS02, HDIS02. 
For models T4\ TIS04, and HDIS04, in which 256 grid cells 
were used in the azimuthal direction, the number of cells in the 
vertical direction was decreased by a factor of two. This is to 
save computational time. Even with this reduction of N z , model 
TIS04 required more than 1100 hours of CPU time per orbit 
(measured at the initial outer edge of thi disk) on a Pentium 
3.06GHz Xeon chip. This is mostly because of the large num- 
ber of Fourier components used in that case (m„ wx = 128 in 
both HDISC4 and TIS04). 

3. Results 

3.1. Hydrodynamic simulations 
3.1.1. Disk fragmentation 

I first present the results of the fiducial hydrodynamic model 
HDIS03, for which N^ = 128. After about one orbit, a strong 
gravitational instability develops in the disk as a consequence 
of its large mass. The early evolution of the instability is sim- 
ilar to what was obtain ed using an adiabatic equation of state 
jFromang et aljE004bl) . However, the spiral arms that form in 
the present case are stronger because of the weaker pressure 
support. When the gravitational instability saturates, the disk 
fragments into dense clumps. This can be seen in figure ^ 
which shows six snapshots of the logarithm of the density in 
the equatorial plane of the disk at times (measured in the or- 
bital time at the initial outer edge of the disk) 9.0, 9.2, 9.3, 
9.5, 9.9 and 10.3 (from top left to bottom right). Note that the 
results of the simulation have been extended by symmetry to 
cover the range [0, 27r] in this figure, which explains why an 
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t(orbits) 


9.5 


9.9-1 


9.9-2 


10.3 




n (so 






07 


pmax 


1UU 








M , ,^ / Mt~ 

IV 1 [ Q I 1V1Y 


3.7 x 10~ 3 


2.8 x 10~ 3 


4.4 x 10~ 3 


2.4 x 10~ 3 




O.D X 1U 


J.o X 1U 


7 A w i n-3 


4.1 X 1U 


MjjMj 


6.5 x lfr 4 


4.3 x 1(T 4 


4.4 x 10~ 4 


6.2 x lO -5 


' clump 


2 x ltr 2 




2 x 10~ 2 


2 x 10- 2 


r H 


8.3 x 10~ 2 


5.7 x 1(T 2 


8.3 x 10~ 2 


0.1 


a 


0.33 


0.5 


0.28 


0.11 



Table 2. Properties of the four clumps visible in the bottom row 
of figure ^ ( m °del HDIS03). The first line of the table shows 
the corresponding time at which each of these properties are 
calculated. It serves as a way of labelling the clumps. The fol- 
lowing lines then give the orbital radius a of the clump, its peak 
density p max , its mass M\q (resp Mioo), ie the mass contained 
in the volume where the density is larger than p max /10 (resp 
Pma.i/100), normalized by the total mass of the disk Mr, the 
Jeans mass My, the radius r c \ ump of the clump, its Hill radius 
and finally the ratio a between its thermal and gravitational en- 
ergy. 



3.0X10" 4 : 

8 2.5x10 -4 - 
c 




-5 nxin~ 5 F - 

7 8 9 10 11 

time 



Fig. 3. Time evolution of the volume averaged gravitational 
stress tensor for model HDIS02 (dotted line), HDIS03 (solid 
line) and HDIS04 (dashed line). The peak value of (TZ av )(t) 
is comparable in all models even though their resolution in the 
azimuthal direction is different. 



even numbers of clumps is always seen. Figure [0 clearly illus- 
trates the fact that the state of the disk is changing very rapidly. 
The clumps form, their density increases, reaches a maximum 
and starts to decline. They eventually completely dissolve in 
the background flow in less than the dynamical timescale. They 
barely complete more than a quarter of their orbit before disap- 
pearing. 

Table [2] presents a list of the properties of the fragments 
at times 9.5, 9.9 (2 clumps exist in the disk at that time and 
are labelled 9.9-1 and 9.9-2) and 10.3, when 3D outputs of 
all the variables were saved. The first and second lines show 
the radius a at which each clump has formed and its maximum 
density p max . The mass Mio and Mioo, normalized by the total 
disk mass M7- are given by the third and forth lines. They re- 
spectively correspond to the mass of the gas contained in the 
volume where p > p„ mx /10 and p > p max /lOO. This should 
be compared with the local Jeans mass, given by the fifth line. 
Line 6 and 7 respectively give the radius of the clump r c [ ump 
and its Hill radius r#, defined by 



where M c — 1 is the central point mass. Finally, the ratio a 
between the thermal and the gravitational energy of the clump 
is given by the last line. 

For each cases, the maximum density is above 100 (the 
maximum density in the initial disk model is 1). The mass of 
the clumps normalized by the total disk mass is a few times 
10~ 3 . Since the latter is of the order of the central point mass, 
it means that the fragments have masses of the order of the 
mass of Jupiter (for a solar-type central mass). The following 
lines of Table [2] give some insights into the stability proper- 
ties of the clumps: their mass is about one order of magnitude 
larger than the local Jeans mass. Their radius is smaller than the 



Hill radius: they are not destroyed by tidal effects. Finally, a is 
smaller than one, showing that they are gravitationally bound. 

All these parameters, taken together, indicates that the 
fragments should not disappear. But still, they dissolved very 
rapidly in the surrounding disk. This is detailed in figure|2 for 
the clump that formed around t - 10.3 (bottom right panel 
of figure Q, density contours in the vicinity of the fragment 
are represented. The maximum of the density in the equatorial 
plane of the disk decreases from 93 to 3.6 between t = 10.21 
and t = 10.37, during which the fragment performed less than 
one fifth of its orbit. The reason for this disappearance, most 
likely numerical, is discussed in section l4~T1 

3.1.2. Numerical effects 

Given the difficulties in accurately describing fragmentation 
in protoplanetary disks, it is important to investigate the nu- 
merical artifacts that could influence the outcome of the sim- 
ulation presented above. The most important of such effects 
is obviously numerical resolution, especially in the azimuthal 
direction. In model HDIS02, N$ equals 64, while in model 
HDIS04, = 256. 

The qualitative evolution of model HDIS04 is very simi- 
lar to that of model HDIS03: fragments form rapidly after the 
gravitational instability saturates. The maximum density they 
reach in larger (p max ~ 500), but, like model HDIS03, they dis- 
appear very quickly. For the low resolution run HDIS02, how- 
ever, no fragment appears at all in the disk. This trend shows 
how sensitive the issue of resolution is on disk fragmentation. 
Indeed, this qualitative difference arises even though the de- 
velopment of the instability seems to be very similar for all 
models. This can been seen in figure [5] which shows the time 
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Fig. 2. Detailed structure of the clump shown in the bottom right frame of figure^at times t = 10.21, 10.29 and 10.37. p max 
is respectively 93,70 and 3.6 at these times. Six contours levels are represented in each panels. They correspond to the density 
values: logp = logp max - 0.5, logp max - 1, etc... 





Fig. 4. Logarithm of the density in the equatorial plane of a 
turbulent disk (model TIS03) at time t = 6.3, when the density 
reaches its largest value p max = 16. 



Fig. 5. Same as figureg] but for model TIS04 at time t = 6.2. 
The largest value of the density in this case is p max = 550. 



history of t he volume averaged g ravitational stress tensor. It is 
defined as jFromang et all2004bl) : 



1 

4^G 



d<S. ; d<fr, 
~di r ~d0 



drdOdz ■ 



(10) 



The time history of (T^ av ) is shown in figure [3] for model 
HDIS02 (dotted line), HDIS03 {solid line) and HDIS04 
{dashed line). For each cases, it starts to rise after about 8.5 or- 
bits, and reaches similar peak values of 2.5 x 10~ 4 . The curves 
are all very similar, showing that the azimuthal resolution has 
little effect on the strength of the instability itself, despite the 
qualitative differences in the fragmentation of the disk. 



disk fragmentation is observed, while fragments quickly form 
in the high resolution model TIS04. This is illustrated by fig- 
ure 0] and |5] The former shows the logarithm of the density 
in the midplane of the disk for model TIS03 after 6.2 orbits, 
which corresponds to the time when it reaches its largest den- 
sity, p max ~ 16. Note how different the disk structure looks 
compared to its hydrodynamic equivalent HDIS03 (figure^. 
On the contrary, figure [5] clearly shows that fragments formed 
in model TIS04, despite the presence of MHD turbulence. The 
properties and time history of the clumps are very similar to 
those observed in the hydrodynamic models described above. 
The likely origin of the difference between models TIS03 and 
TIS04 is discussed in the following section. 



3.2. MHD simulations 

In this section, I compare the results of model TIS03 and 
TIS04 with their hydrodynamic counterpart HDIS03 and 
HDIS04. The results of model TIS02, for which = 64, are 
very similar to those of model TIS03 and will not be described 
further. 

In both TIS03 and TIS04, a gravitational instability 
quickly develops on top of MHD turbulence. However, the 
qualitative evolution of the instability is dramatically differ- 
ent in both models. In the lower resolution case TIS03, no 



4. Discussion 

4.1. The Jeans length 

A critical parameter to take into account in the simulations pre- 
sented above is the Jeans length, defined by: 



Gp 



(11) 



Truelove et al 1 dl997h showed that spurious fragmentation can 
occur if the Jeans length is not well resolved at all times during 
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1 ,000 f 



0.100 - 




0.001 I , , , I , , , I , , , I , , , I , , , I , , , 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 

r 

Fig. 6. Jeans length vs radius for model HDIOS3 at time t — 
{dotted-dashed line) and t = 9.9 (solid line). The dashed line 
shows the maximum of the radial and azimuthal spacing of the 
grid as a function of radius. The dotted line plots four times this 
value. 



3.0X10" 4 T 

8 2.5x10 -4 - 
c 




o 

-5 oxi rr 5 F - 

5 6 7 8 9 10 11 

time 



Fig. 8. Time history of the volume averaged gravitational stress 
tensor for models HDIS03 [solid line), TIS03 (dotted line) and 
TIS04 [dashed line). The earlier growth of the gravitational 
instability in the MHD cases is triggered by the large density 
perturbations due to MHD turbulence. 




0.2 0.4 0.6 0.S 1.0 



Fig. 7. Same as figure[6] but for models TIS03 (dotted-dashed 
line) at time 6.3 and TIS04 at time 6.2. The dotted and dashed 
lines have the same meaning as for figure [6] (the upper curves 
show the case of model TIS03, for which = 128, while the 
lower curves correspond to model TIS04). 



the disk evolution and that the largest grid cell on the com- 
putational domain should be about 4 times larger than Aj. In 
isothermal simulations, it is very difficult to resolve the Jeans 
length during the disk fragmentation: since the sound speed is 
constant, Aj increases as (1/p) 1 ^ 2 . It is therefore important to 
follow its behaviour. At each radii, one can calculate the mini- 
mum value of Aj in the equatorial plane. Its radial profile is rep- 
resented in figure|5|for model HDIS03 at times t = (dotted- 
dashed line) and t = 9.9 (solid line). The dashed line represents 
the largest grid size as a function of r and the dotted line shows 
4 times this value. At t — 0, it is obvious that the Jeans length 
is well resolved and satisfies the Truelove criterion. However, 
at t = 9.9, two fragments formed in the disk and the resolution 
of the simulation becomes coarser than the local Jeans length. 
According to the Truelove criterion, the evolution of the frag- 



ments after this point is not meaningful and should not be re- 
garded as being physical. On top of this, the size of the forming 
clumps themselves at this stage of the simulation is of the order 
of a few grid cells. The gravitational potential they create (and 
which governs their subsequent collapse) is not accurately eval- 
uated at this scale because of the finite resolution of the grid. 
Their later evolution in the simulation is not well described. 
Both of these arguments show that these calculations can only 
be used as a diagnostic for disk fragmentation. 

As a conclusion, the hydrodynamic simulations tell us that 
massive quiescent disk should fragment in the conditions stud- 
ied here. However, they cannot tell us anything about their long 
term evolution, even though their stability properties, presented 
above, suggest that they may survive. 

Figure is similar to figure [6] it compares the grid reso- 
lution with the local Jeans length for the MHD turbulent disk 
models TIS03 (dotted-dashed line) and TIS04 (solid line) re- 
spectively at times t = 6.2 and t = 6.3 (ie at the same time as 
the snapshots shown in figure |4] and [5}. In the case of model 
TIS03, the minimum value of the local Jeans length (which 
decreases at the location of the spiral arms), is roughly three 
times the resolution of the simulation. The Truelove criterion 
is only marginally violated in that case. However, in similar 
conditions, clumps are seen to form in model HDIS03, which 
shows that MHD turbulence clearly has an effect in that case. 
FigureQalso shows that the case of model TIS04 is very sim- 
ilar to the hydrodynamic models for which fragmentation is 
found: the local Jeans length decreases as the density of the 
fragments increases, until the former becomes smaller than the 
resolution. The subsequent evolution of the clump is not accu- 
rately described by the simulations. 
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4.2. The effect of MHD turbulence 

The results of model TIS03 and TIS04, together with those of 
their hydrodynamic counterpart HDIS03 and HDIS04 suggest 
that MHD turbulence affects the outcome of gravitational insta- 
bilities in disks. Although the precise interplay between these 
processes is difficult to understand given the complexity of the 
problem and its sensitivity to numerical resolution, some indi- 
cations are provided by a detailed analysis of the simulations. 
Figure[S]shows the time history of the volume averaged gravita- 
tional stress tensor (Tf?™) (defined in section l3~.1.2> for models 
HDIS03 (solid line), TIS03 (dotted line) and TIS04 (dashed 
line). The maximum value reached by (T*™ v ) is smaller when 
the disk is turbulent, suggesting that the strength of the gravita- 
tional instability itself is reduced by the presence of the turbu- 
lence. It may explain why no fragments form in model TIS03, 
as the gravitational stress is decreased by roughly a factor of 
two compared to model HDIS03 (Note that, even though the 
evolution of the stress is shown for a rather short period of 
time, the difference between the various runs is meaningful be- 
cause this is precisely during that time that fragments forms 
or not). This effect of MHD turbulence on the strength of the 
gravitational instability was al ready noted in s imulations using 
an adiabatic equation of state IFromang et alj|2004bi) . In adia- 
batic simulations, this decrease of the gravitational stress did 
not affect the qualitative outcome of the gravitational instabil- 
ity. However, in model TIS03 it prevents the density to become 
large enough in the spiral arm to trigger their collapse. 

Figure [8] also shows that (Tf™ v ) is larger in model TIS04 
(in which clumps form) than in TIS03. This is not due to the 
gravitational instability alone, since no such dependency is ob- 
served in the hydrodynamic models HDIS02, HDIS03 and 
HDIS04. Nor is it due to the MHD turbulence alone, since its 
properties are similar in model TIS03 and TIS04. However, 
as is explained below, the disk fragmentation (as well as the 
increase of the gravitational stress tensor) observed in model 
TIS04 could be due to the small scale angular transport prop- 
erties of the MRI. In cases, such as the present one, when the 
magnetic field is mostly toroidal, the most unstable mode of 
the MRI has an azimuthal waven umber m cr „ that satisfies the 
relation (Balb us & Hawlevll998l) 



To summarize, these results suggests that, on top of the 
Jeans length, it is important in calculations of the present type 
to well resolve the length Ab defined by 



R 

which gives 




(12) 



(13) 



At the position where the fragment forms in model TIS04 
(see figure the properties of the disk model indicate that 
R£l/c s ~ 16, while the turbulence is characterized by f3 ~ 10. 
These values gives m cril ~ 35. Such a mode is barely resolved 
in model TIS03 (only 7 grid cells per wavelength) but much 
better described in model TIS04 which has twice the resolu- 
tion. In the latter case, it grows on an orbital timescale and helps 
to remove the angular momentum of the collapsing body. In the 
former case, the growth rate of this small scale mode is lower 
and it is likely it will be less important in the collapse process. 



A B = v a P 



(14) 



where P is the orbital period in the disk. 

Final ly, it is interestin g to compare the present results to 
those of lKim et alJ ( 120031) who performed local calculations of 
giant molecular cloud formation in the presence of MHD tur- 
bulence. They found that density fluctuations triggered by the 
MRI can be swing-amplified and lead to disk fragmentation. 
This is not observed here. The reason for this d i fferenc e lie 
in the initial setup of both calculations: iKim et alJ J2OO3I start 
their simulations with vertical fields with a non-zero net flux 
that produces a stronger turbulence than in the present work. 
Indeed, the ratio of the total stress to the pressure they ob- 
tain is ~ 0.2, while it is an order of magnitude smaller here 
dFromang et all2004 b). The larger density fluctuations that re- 
sult are more likely to be gravitationnally bound in their case 
than in the present one. 

5. Conclusion 

This paper presents a collection of hydrodynamic and magne- 
tohydrodynamic 3D numerical simulations of the evolution of 
massive isothermal disks. 

In hydrodynamic simulations, the disk is found to fragment 
into high density clumps when the number of grid cells N$ 
is larger than 128 in the azimuthal direction. The fragments, 
whose mass is comparable to the mass of Jupiter, disappear 
very quickly in the disk. However, due to the limited resolu- 
tion of the simulations, the local Jeans length is no longer re- 
solved when their density has increased by about two orders 
of magnitudes. Moreover, at these small scales, the collapse 
of the clumps themselves is not accurately described because 
the gravitational potential is smoothed at the grid scale. This 
is probably what causes their quick disappearance. At the very 
least, the long term behavior of these fragments cannot be stud- 
ied with these simulations, although various criteria suggest 
that they could survive for an extended period of time. 

When a weak magnetic field is present initially in the disk, 
it becomes turbulent because of the nonlinear development of 
the MRI. This turbulent nature of the flow has an effect on the 
evolution of the gravitational instability: no fragment is found 
when Afy = 128, in contrast with the hydrodynamic model hav- 
ing the same resolution. However, fragments appear when the 
azimuthal resolution is increased to Afy = 256. An analysis of 
these results suggest that they are a consequence of two com- 
peting effects: on the one hand, MHD turbulence tends to de- 
crease the strength of the gravi tational instability, as was al- 
ready noted in a previous study jFromang et al.ll2004bl) : this is 
sufficient to suppress gravitational fragmentation in the lower 
resolution case. On the other hand, when the most unstable 
modes of the MRI are resolved, small scale angular momentum 
transport appears to help the formation of bound fragments. 
This interplay between disk fragmentation and radial transfer 
of angular momentum within an overdense region was already 
suggested bv lKim et alJ |2003). 
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Despite the difficulties related to the limited numerical res- 
olutions of the simulations, the results presented in this paper 
indicate that MHD turbulence is an important component in 
the process of massive disk fragmentation. Future studies are 
needed to quantify the importance of the turbulence and to in- 
vestigate whether the fragments that form are truly long lived or 
quickly dissolve in the background flow. To cope with the large 
resolutions needed to adress this problem, adaptative mesh re- 
finement techniques, possibly in cylindrical coordinates, may 
prove to be the method of choice. 
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